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We study the angular resolution of the gravitational wave detector LISA and show that numerical 
relativity can drastically improve the accuracy of position location for coalescing Super Massive 
Black Hole (SMBH) binaries. For systems with total redshifted mass above 10 7 Mq , LISA will mainly 
see the merger and ring-down of the gravitational wave (GW) signal, which can now be computed 
numerically using the full Einstein equations. Using numerical waveforms that also include about 
ten GW cycles of inspiral, we improve inspiral-only position estimates by an order of magnitude. 
We show that LISA localizes half of all such systems at z = 1 to better than 3 arcminutes and the 
best 20% to within one arcminute. This will give excellent prospects for identifying the host galaxy. 
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Introduction. — Astronomical observation provides 
very strong evidence for the existence of a "dark" com- 
pact massive (10 6 — 10 9 Mq) object in the core of every 
galaxy for which the central parsec region can be resolved 
PP. These objects are believed to be supermassive black 
holes (SMBH) , and are of great interest to researchers in 
fundamental physics, astrophysics and cosmology. Their 
formation and the observed correlation between SMBH 
mass and galaxy morphology (see [5] for an overview) are 
still open questions. SMBHs probably arise at least in 
part through the mergers of smaller-mass BHs |3J . These 
mergers and the mergers of SMBH binaries following col- 
lisions of galaxies constitute some of the most powerful 
sources of gravitational waves (GW) predicted by current 
models. We will be able to detect such events through- 
out the Universe with LISA pi], a proposed space-borne 
gravitational-wave observatory, scheduled for launch in 
2018+ and designed to be sensitive to GW signals in the 
range 1(T 4 - 0.1 Hz. 

Realizing the full scientific benefit of the LISA mission 
will require accurate estimates of the binary's parame- 
ters. Precise measurements of the masses, spins and dis- 
tance will allow us to probe models of SMBH formation. 
Accurate localization of the source in the sky is crucial 
to relate gravitational-wave and electromagnetic obser- 
vations of the coalescence event, and hopefully will allow 
identification of the host galaxy. Optical observations 
are required to measure the redshift to the object, while 
gravitational- wave observations yield precise calibration- 
free estimates of the distance. Such LISA events with 
optical counterparts will determine the redshift-distance 
relationship, which in turn will allow us to map the ge- 
ometry of the Universe and measure the amount of dark 
energy. 

Recently several groups [SJ [BJ [5] have evaluated the 
accuracy of parameter estimation using the inspiral part 
of the GW signal. It was shown that the spin and higher 
orbital harmonics are necessary to de-correlate the pa- 
rameters and therefore improve the parameter estima- 



tion. 

In this Letter we assess the angular resolution of 
LISA for SMBH binaries with (red-shifted) masses above 
10 7 Mq. We expect several such merger signals per year 
at relatively close distance [3]. For those heavy systems 
the inspiral signal may be smaller or at least not much 
larger than the instrumental noise, and the signal will 
be dominated by the merger and ring-down. We use 
numerical relativity to compute a waveform which con- 
tains about ten GW cycles, plus merger and ring-down. 
We fix the redshifted masses of two non-spinning BHs to 
4.44 x 10 6 M Q and 8.88 x 10 6 Mq and vary the "extrinsic" 
parameters: sky ecliptic coordinates 9s, </>s, inclination 
i of the orbital angular momentum to the line of sight, 
polarization angle vb, fiducial arrival time To which we fix 
to be the time when the binary separation equals 10 M, 
orbital phase 4>q at T , and luminosity distance D^. 

Data analysis for LISA is based on time-delay in- 
terferometry (TDI, see [9] for an overview). We con- 
vert the strain polarizations h + and h x given in the 
source frame to first generation unequal-arm Michel- 
son streams X, Y, Z [TTJ] and use two combinations 
A = (2X - Y - Z)/3, E = (Z - Y)/V3 with uncorre- 
lated noise. Due to technical difficulties explained below, 
we do not take into account the third independent combi- 
nation, which has poor sensitivity to GW at low frequen- 
cies, and which adds only a few percent to the combined 
signal-to-noise ratio (SNR). By fixing masses we under- 
estimate the errors, but at the same time not including 
the third TDI combination overestimates the error boxes. 
We mainly concentrate here on the estimation of the lo- 
calization of the source and usually 8 the directional 
angles very weakly correlate with masses. Since we use 
the merger for localizing the hosting galaxy, we cannot 
produce an early warning for the merger itself; however, 
we can identify the hosting galaxy by the afterglow [TT]. 

We conduct Monte Carlo simulations by randomly 
choosing 600 points in the extrinsic parameter space 
(with fixed masses, and choosing an example distance 



2 



of z = 1 or Dl ~ 6.4 Gpc), and estimate the errors in 
the parameter error box with two different methods: the 
first is based on computing the variance-covariance ma- 
trix, the second is a Bayesian method based on the eval- 
uation of the marginalized posterior distribution func- 
tion using a Metropolis-Hastings Markov chain (MHMC) 
[T2l [T3] E] • We have shown that both methods give com- 
parable results. For 50% of the randomly chosen parame- 
ters we can localize the source down to a box 3 arcminutes 
on a side, and the best 20% of the events are localized 
to better than one arcminute. The relative error in the 
distance is less than 1%, the largest error in the distance 
being due to weak lensing [15] . We also estimate robust- 
ness of our results with respect to errors in the numerical 
waveform. By conducting another Monte Carlo simula- 
tion, we have found that the errors in the sky locations 
are good to within factor two at worst (but usually much 
better than that). 

Numerical relativity waveforms. — We use numerical- 
relativity waveforms from simulations of non-spinning 
black holes at mass ratio 1:2, which have initially been 
presented in |16j . and are discussed in more detail in |17j . 
where the quadrupole spherical harmonic mode has been 
compared with effective-one-body waveforms, and excel- 
lent agreement in the phase evolution has been found. 
Here we now include results for higher spherical-harmonic 
contributions. Our simulations follow the now stan- 
dard moving-puncture approach [T5J , using the BAM 
code [20, 21 to numerically solve the Einstein equations. 
Black hole initial data are modeled as conformally flat 
puncture data [22, 23 . The initial data parameters then 
reduce to the black-hole masses, momenta and separa- 
tion. We choose a separation of 10 M (we refer to the 
total initial black-hole mass as M). The momenta are 
specified to give quasicircular inspiral with minimal ec- 
centricity of e w 0.003 [24] . 

The gravitational-wave signal is extracted at five 
surfaces of constant radial coordinate, R ex = 
40,50,60,80,90M by means of the Newman-Penrose 
Weyl tensor component ^4, as described in [5D]. At ev- 
ery extraction radius the gravitational wave strain is ob- 
tained from ^4 by double time integration as described 
in |16j . The analysis carried out in this paper will use, 
as approximate asymptotic amplitude, the curvature per- 
turbation extracted at radius 90M, at our highest numer- 
ical resolution, as has been done in [T7]. For the modes 
with I > 2 we filter out low frequencies and frequencies 
higher than the Nyquist frequency to avoid the patholo- 
gies discussed in Sec. II. A. of [55] ■ For these simulations, 
we find that the finite extraction radii dominate the er- 
ror, and the amplitude error is below 5% prior to merger, 
and the accumulated phase error is below 0.25 radians for 
the 700 Af up to Mu = 0.1. The fall-off in the amplitude 
error with respect to radiation extraction radius is not 
clean around merger time, preventing us from performing 
an accurate extrapolation to infinity. As such, we would 



conservatively give an uncertainty estimate of 10% of the 
amplitude at merger and later. We estimate the relative 
error in the amplitude between different modes (which is 
what dominates the results presented here) as below 4%. 

Parameter estimation. — We have used two methods 
to estimate accuracy in measuring the parameters of the 
GW signal. The first method is based on computing the 
variance-covariance matrix. This method is widely used 
and well described in numerous publications (see for ex- 
ample [15, 26, 27] [55]). It is based on inverting the Fisher 
matrix T, which is the matrix of the inner products be- 
tween derivatives of the signal with respect to the param- 
eters h t i = dh/dXi: 

r°° hi(f)h*Jf) 

r« = 4« J df ' g n(/ y , (o-i) 

where S n (f) is the one-sided noise power spectral den- 
sity. We use two TDI measurements and the combined 
Fisher matrix is a sum of Fisher matrices for A and E: 
Tij = Tfj + Tfj . The presence of the noise might cause 
a deviation of the recovered parameters of the GW sig- 
nal from the true values. The diagonal elements of the 
variance-covariance matrix are maximum likelihood es- 
timators of the variance of parameters around the true 
value in the case of a large SNR (which is always the 
case here). We have computed derivatives numerically 
using forward differencing and checked the robustness 
with respect to the step using several randomly chosen 
points in the parameter space. We have tried to choose 
a parametrization of the signal which would reduce the 
dynamical range of the eigenvalues of the Fisher matrix, 
however there are still six or seven orders of magnitude 
between the smallest and largest eigenvalues. 

As a second way to evaluate the parameter estima- 
tion error we use a Metropolis-Hastings Markov chain 
(MHMC) [12] to estimate the posterior probability func- 
tion 

p(X\s) cx tt(A)Z(s|A), (0.2) 

where 7r(A) is the prior distribution of parameters and 
L(s\X) is the likelihood function. The MHMC uses a pro- 
posal distribution g(A(j)|A(j\) and Metropolis-Hastings 
ratio 

H = P(A(j)|a)g(A (3 -)|A w ) 
p(A (i )|s)g(A w |A (j )) ' 

which gives a probability a(A(j)|A(j)) = mm(l, H) of ac- 
cepting the jump from the point A/j) to X(j)- MHMC 
gives the best result (better convergence) if the pro- 
posal distribution matches the shape of the target dis- 
tribution. For the proposal distribution we take nor- 
mal jumps in the eigen-directions of the Fisher ma- 
trix |14j . yielding acceptance rates > 30%. We have 
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generated simulated data using the lisatools software 
[http://code.google.eom/p/risatools7 , consisting of in- 
strumental noise and a simulated reduced Galaxy con- 
fusion noise with more than 5 x 10 7 chirping binaries 
[2UI [3D]. We have generated signals for 600 randomly 
chosen parameters (600 data sets), and performed 600 
mappings of posterior distribution functions using 3 x 10 5 
long chains. We have also used mild simulated anneal- 
ing for the first two thousands steps, which helps the 
chain to migrate from the point with true parameters to 
the point with better likelihood. Note that due to the 
presence of the Galaxy the noise is strictly speaking not 
Gaussian. Because the jumps are rather small, we have 
made the assumption that the Fisher matrix does not 
change notably within the jumping region and therefore 
— an d the Metropolis-Hastings 

ratio is reduced to Metropolis. This assumption results 
in a small bias in the estimation of the variance, but it 
significantly speeds up the computations. 

The results of both methods are presented in the top 
two plots of Figure [T] as cumulative histograms for the 
600 realizations. Both methods give comparable estima- 
tion of the error box for the sky location (0s,<j>s)- One 
can see that 50% of the trials give an error box smaller 
than 3 arcminutes for a source located at z = 1. 
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FIG. 1: Top plots: cumulative histograms for the standard 
deviation estimators from 600 randomly chosen points in the 
parameter space. The black (lower) curve corresponds to the 
variance a evaluated by computing the variance-covariance 
matrix and the red (upper) curve is the Bayesian estima- 
tor obtained using the MHMC. Lower plots: deviation of the 
mean value of the chain from the true parameter of the sig- 
nal. The plotted results correspond to a distance of 6.4 Gpc 
(z » 1 and total mass of « 0.7 x 10 7 M©), 

The duration of the signal from a separation of 10M 
to the merger is 26 hours, so the Doppler modulation is 
small, however for our choice of parameters in the Monte 
Carlo simulation the SNR varies between 900 and 9000, 



and these large values drastically improve parameter es- 
timation. The other crucial ingredient for the excellent 
angular resolution is the presence of higher orbital har- 
monics coming from the higher multipoles of the source. 
We find that the use of only the dominant (I — \m\ = 2) 
mode worsens our results by up to a factor of ten. Dif- 
ferent harmonics have different angular emission patterns 
and different dependencies on the inclination. These help 
to de-correlate the parameters, and together with the 
high SNR turn a seemingly small effect (higher modes 
are much lower in amplitude) into a crucial contribution 
for parameter estimation. For our analysis we have used 
Z = 2,3,4 and m = —I, . . . I (except m = 0). We also find 
that the arrival time could be measured with an accuracy 
of less than a second and the luminosity distance with an 
accuracy less than 1.5%. However, the dominant error in 
estimating the luminosity distance is due to weak lensing 
[l"5] . Weak lensing gets stronger with distance [31] and 
we expect to have SMBH binaries at distances z < 5. 

As we have mentioned, the Markov chain is usually mi- 
grating to the point with higher likelihood and this shift 
could be as large as the variance itself. We have found 
that the secondary maxima in the likelihood could be 
very close to the primary maximum both in amplitude 
and distance in parameter space. In our simulations we 
observed that due to the presence of the instrumental and 
Galactic confusion noise the secondary maxima could be- 
come stronger than the primary (by less than 0.01%) and 
the secondary could be located about a a away, where a 
is the standard deviation of the parameter assuming no 
secondary maxima. In the lower two plots of Figure [T] we 
show the histogram of the deviation mean value of the 
chain from the true sky location. Further work is needed 
to explore ways of reducing these ambiguity errors. 

We have used only two TDI streams in the above 
analysis, because we could not compute the third inde- 
pendent data set constructed out of X, Y, Z, which is 
T cx X + Y + Z, with the required accuracy: the sig- 
nal cancellation at low frequencies was not perfect due 
to inaccuracy (less than a few percent) in the evalua- 
tion of X, Y, Z. In principle we could generate A, E, T 
out of the six-pulse TDI combinations as in [32] but we 
did not have simulated a Galactic background for those 
combinations, so we have decided to drop the third effec- 
tive detector for the present considerations. We however 
checked that using T (as in [35]) would improve our SNR 
by a few percent, and, more importantly, could also im- 
prove parameter estimation: having three independent 
measurements should help to triangulate the source. 

Finally, we have checked the robustness of the vari- 
ance estimation with respect to possible errors in the 
numerical waveform. Since the results depend on the rel- 
ative strengths of the different harmonics, we varied the 
I = m = 2 mode by ±5% with respect to the higher har- 
monics. We have estimated the variance by computing 
the variance-covariance matrix and comparing it to the 
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original one, and find that a changes by at most a factor 
of two (much less in the majority of the cases < 15%). We 
have also noticed that changes in the waveform resulting 
from enhancing the higher orbital modes improve results 
(makes the variance smaller) which is in agreement with 
our explanation of the angular resolution. 

Summary. — We have presented a first study of how 
results from numerical relativity can improve parame- 
ter estimation for SMBH binary observations with LISA. 
We have chosen to first examine the issue of angular 
resolution of LISA - which is crucial to identify elec- 
tromagnetic counterparts to gravitational wave observa- 
tions. Looking at the merger signal for non-spinning 
SMBH binaries with redshifted total mass 1.4 x 1O 7 M 
and a mass-ratio of 2, Placing our source at a luminos- 
ity distance of 6.4 Gpc (z « 1, corresponding to a total 
mass of « 0.7 x 10 7 A/ Q ), we have found an excellent 
sky resolution: for 50% of our randomly chosen events 
we can localize the source down to 3 arcminutes, which 
roughly corresponds to an order of magnitude improve- 
ment over estimates using only the inspiral phase of the 
GW signal. We have also shown that we obtain con- 
sistent error estimates when comparing two independent 
methods based on the variance-covariance matrix and a 
Metropolis-Hastings Markov chain. 

Our result calls for further work along several direc- 
tions. First, it is necessary to cover the parameter space 
of numerical simulations: larger mass ratios will show 
two effects: a decrease in SNR, but an increase of the 
contributions of higher modes, which have proven crucial 
for parameter estimation. The inclusion of spins will lead 
to a much larger parameter space, but also to much more 
interesting waveform structures. For lower masses, it will 
be important to accurately match numerical waveforms 
to post-Newtonian results in order to cover the whole 
LISA band, and to understand the systematical errors in 
the matching procedure. The extension of existing phe- 
nomenological waveforms (e.g. along the lines of [T6|l33j) 
to spinning binaries and non-dominant modes will allow 
for systematic parameter estimation studies in the regime 
where solving the full Einstein equations is necessary to 
obtain accurate gravitational wave signals, and will ex- 
tend the useful sensitivity range of LISA to substantially 
higher masses, so that LISA can begin to explore the pro- 
cess whereby black holes grow from the size we see in our 
galaxy to the sizes that are needed in quasars. 
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